Motivation of the study

This study aims to calibrate Leaf Area Density (LAD), which is defined as the area of leaves per unit of tree crown volume (in \(m^2.m^{-3}\)). This parameter forms the basis of numerous process-based forest models, such as physiological models that consider photosynthesis and transpiration (e.g. SurEau by Cochard et al., 2021), and forest gap models that estimate light interception (e.g. ForCEEPS by Morin et al., 2021, and Samsara2 by Courbaud et al., 2015).

Although the leaf area density can be directly estimated from a sample of trees (Stadt & Lieffers, 2000), this method is laborious and difficult to replicate for every tree in a stand. Therefore, indirect methods can be employed, such as model inversion using a hemispherical photograph (Courbaud et al., 2003) or the more recent LiDAR methodology (Wei et al., 2020). However, there is a significant lack of research on LAD in the scientific literature, as there are no large-scale databases containing values for a wide range of European species. Furthermore, the few reported species-specific values vary considerably between studies (Ligot et al., 2014).

The aim of this study is to estimate leaf area density (LAD) for nine different species/groups: Fagus sylvatica, Picea abies, Quercus sp., Carpinus betulus, Pseudotsuga menziesii, Pinus sylvestris, Abies alba, Larix decidua and ‘other species’. Additionally, Nock et al. (2008) demonstrated intraspecific variation in LAD, with intracrown leaf area decreasing by up to 40% with increasing tree age. This led us to consider the effect of diameter as well. Finally, as crown dimensions are very sensitive to local competition (Touzot et al., 2025), we hypothesise that LAD may also vary in response to competition.

Methodology

In this study, we used a model inversion methodology, using the ray-tracing model SamsaraLight (Courbaud et al. 2003). We want to estimate the SamsaraLight tree-level parameter LAD using measurements of transmitted light \(PACL_{obs}\) (for percentage of above light canopy), based on hemispherical photographs, within 45 plots of different composition and structures. For each plot, trees are inventoried with their crown dimensions, allowing to derive an estimated transmitted light \(PACL_{est}\) on virtual sensors with the ray-tracing model SamsaraLight. By doing so, for each virtual sensor of each plot, we have the observed transmitted light \(PACL_{obs}\), and the estimated transmitted light given a tree-level value of LAD \(PACL_{est}\).

The aim of the calibration is to find a tree-level LAD model that minimize the error between \(PACL_{obs}\) and \(PACL_{est}\). To do so, we defined the LAD for a tree \(i\), of species \(s\) in a plot \(p\) as \(LAD_{i,s,p} = a_{0,s} + a_{0,p} + a_{1,s}.DBH_i + a_{2,s}.BAH_i + \epsilon\), where \(a_{0,s}\) is a species-specific intercept, \(a_{0,p}\) is a random hierarchical site effect, nested within an origin effect (4 origins of the plots within UCLouvain, INRAe, Cloture and IRRES), \(DBH_i\) is the diameter at breast height of tree \(i\) , with the associated estimated species-specific coefficient \(a_{1,s}\), \(BAH_i\) the competition index of the tree \(i\) (i.e. the sum of basal area of trees higher than the tree \(i\)), with the associated estimated species-specific coefficient \(a_{2,s}\), and \(\epsilon\) the normal error term. We used a Bayesian approach to estimate the posterior distribution of parameters that minimize the log-likelihood between observed and simulated PACL. As computation of \(PACL_{est}\) needs an external call of SamsaraLight model, we used the R package ‘BayesianTools’ to run a MCMC sampling of parameters with the “DREAMzs” sampling algorithm.

For this study, we used the SamsaraLight model with the R package SamsaraLight (Beauchamp et al., work in progress, intended for publication in the Journal of Open Source Software). The SamsaraLight model was initially developed in Java and included in the Capsis simulation platform (SamsaraLightLoader by Gauthier Ligot). It could be called from R software via the “RCapsis” R package (Fortin et al., 2021). To speed up calculation times and make SamsaraLight easier to use in an R environment, we created the “SamsaraLight” R package. The source code is written in C++ to enable fast calculations and communication with the R environment is handled by the Rcpp package (Eddelbuettel & Balamuta, 2018). The “SamsaRaLight” R package reduces calculation time compared to the previous method (the RCapsis R package calling the Java model SamsaRaLightLoader) as it does not require the creation of a Java server nor Java objects. It also facilitates the use of SamsaraLight within R, as it does not require the creation of an inventory file, instead using R objects and functions directly. The package is functional and stored in a GitHub repository, but it is still subject to private restrictions as discussions regarding co-authorship and the intellectual property rights of the source code need to be concluded.

Calibration plots

The 45 calibration plots comes from 4 different sources we wiil call IRRES (9 plots, Gauthier Ligot), CLOTURE (23 plots, Gauthier Ligot), INRAe (10 plots, Philippe Balandier) and UCLouvain (3 plots, Mathieu Jonard).

IRRES

CLOTURE

INRAe

UCLouvain

Preliminary analysis

As a preliminary analysis, we tried to fit a mean LAD value for each sensor of each site. To do so, we run the SamsaraLight model on each stand by fixing the same crown LAD values of each tree, and for each sensor, we computed the residuals between the estimated PACL value on the virtual sensor, and the observed PACL on the field for this sensor. We did so for 500 values of LAD from 0.01 to 5. By doing so, we could find for each sensor the LAD that minimizes the residuals.

Many sensors did not converged (i.e. best fitting LAD was 5 as the residuals function showed an asymptotic form). Indeed, we could increase as much as we want the trees’s LAD, the estimated PACL could not decrease as much as the observed one. PACL is principally linked to the crown geometries, and secondary by the crown LAD. Thus, there are two cases where crown LAD did not influence PACL: (i) for open stands, where the high light on the ground comes principally from unobstructed rays, (ii) for highly crowded stands, where the very low light on the ground comes from the high number of intercepted crowns. Another reason could be linked to the wrong representation of leaves and crown above the sensor, where a simple branch can obstruct the sensor in the field where we considered a simple crown shape that result in incapabilities of representing such a light obstruction with the SamsaraLight model. Consequently, we remove the sensors that did not converged (465/1121 sensors), leading to remove consequently 3 sites where all sensors did not converged: Cloture11, Cloture15, Cloture2.

The mean LAD over all sensors that converged was about 0.944

Bayesian calibration

Clean model: species X DBH with covariance

100,000 iterations and 3 chains

##   id_model n_iterations n_analysis filter_sensors site_rd_effect
## 1        1        2e+05       5000           TRUE           TRUE
## 2        2        2e+05       5000           TRUE           TRUE
## 3        3        2e+05       5000           TRUE           TRUE
##   origin_rd_effect intercept_per_sp intercept_sp_pooling dbh_effect dbh_per_sp
## 1             TRUE            FALSE                FALSE       TRUE      FALSE
## 2             TRUE             TRUE                 TRUE       TRUE       TRUE
## 3             TRUE             TRUE                 TRUE       TRUE       TRUE
##   dbh_sp_pooling consider_covariance compet_effect
## 1          FALSE               FALSE         FALSE
## 2           TRUE               FALSE         FALSE
## 3           TRUE                TRUE         FALSE

Computation time

##   model chain time_hours
## 1     1     1   310.8808
## 2     1     2   317.5449
## 3     1     3   300.6420
## 4     2     1   319.6946
## 5     2     2   317.3686
## 6     2     3   306.4468
## 7     3     1   317.5156
## 8     3     2   301.6323
## 9     3     3   317.2304

Observe convergence

Effect of DBH for each species

## # A tibble: 27 × 6
## # Groups:   model, chain [9]
##    model chain subchain rho_lower rho_median rho_upper
##    <chr> <chr> <chr>        <dbl>      <dbl>     <dbl>
##  1 1     1     1                0          0         0
##  2 1     1     2                0          0         0
##  3 1     1     3                0          0         0
##  4 1     2     1                0          0         0
##  5 1     2     2                0          0         0
##  6 1     2     3                0          0         0
##  7 1     3     1                0          0         0
##  8 1     3     2                0          0         0
##  9 1     3     3                0          0         0
## 10 2     1     1                0          0         0
## # ℹ 17 more rows

Effect of BAxDBH

100,000 iterations and 1 chain

##   id_model n_iterations n_analysis filter_sensors site_rd_effect
## 1        1        1e+05       5000           TRUE           TRUE
##   origin_rd_effect batot_in_site_effect intercept_per_sp intercept_sp_pooling
## 1             TRUE                 TRUE            FALSE                FALSE
##   dbh_effect dbh_interaction_batot dbh_per_sp dbh_sp_pooling
## 1       TRUE                  TRUE      FALSE          FALSE
##   consider_covariance compet_effect
## 1               FALSE         FALSE

Computation time

##   id_model time_hours
## 1        1   54.08484

Observe convergence

Effect of DBH X batot on a mean species

Output analysis

Observe initial data

## 
## Call:
## lm(formula = batot_m2ha ~ order, data = data_output_tree)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.733  -7.807   0.806   5.892  34.167 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      21.4670     0.1306  164.38   <2e-16 ***
## ordergymnosperm  11.6324     0.1889   61.59   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.527 on 8167 degrees of freedom
## Multiple R-squared:  0.3172, Adjusted R-squared:  0.3171 
## F-statistic:  3794 on 1 and 8167 DF,  p-value: < 2.2e-16

Stand variables

Crown dimensions

Tree-level consequences on tree leaf area

Stand-level consequences on leaf area index LAI

Consequences on light

Tree-level

Stand-level